library(survey)
library(foreign)
library(ggplot2)
library(ggrepel)
library(cowplot)
library(dplyr)

## see PCM_cleaning.R for data cleaning
load("~/Dropbox/Public_Conf_Mil/Data_and_code/wave2_data/clean_dataw2.RData")

## Create weighted survey design object
w2_design <-
  svydesign(
    id = ~ 1,
    weights = ~ weight2,
    data = df
  )

#### FIGURE 2.8 - DEMOGRAPHIC FEATURES ####

dm_prty <- as.matrix(svyby(~ Q11_mr, ~ party, w2_design, svymean, na.rm = TRUE))
dm_prty <- cbind(c(rep(15,3)), dm_prty)
dm_prty[,2] <- c("Democrat", "Independent", "Republican")

dm_ideo <- as.matrix(svyby(~ Q11_mr, ~ ideo3, w2_design, svymean, na.rm = TRUE))
dm_ideo <- cbind(c(rep(14,3)), dm_ideo)
dm_ideo[,2] <- c("Liberal", "Moderate", "Conservative")

dm_gend <- as.matrix(svyby(~ Q11_mr, ~ male, w2_design, svymean, na.rm = TRUE))
dm_gend <- cbind(c(rep(13,2)), dm_gend)
dm_gend[,2] <- c("Female", "Male")

dm_race <- as.matrix(svyby(~ Q11_mr, ~ RACETHNICITY, w2_design, svymean, na.rm = TRUE))
dm_race <- cbind(c(rep(8,4)), dm_race[c(1, 2, 4, 6),])
dm_race[,2] <- c("White", "Black", "Hispanic", "Asian")

dm_relg <- as.matrix(svyby(~ Q11_mr, ~ religion, w2_design, svymean, na.rm = TRUE))
dm_relg <- cbind(c(rep(9,4)), dm_relg[2:5,])
dm_relg[,2] <- c("Christian", "Catholic", "Other", "None")

dm_mvet <- as.matrix(svyby(~ Q11_mr, ~ vetstatus, w2_design, svymean, na.rm = TRUE))
dm_mvet <- cbind(c(rep(12,3)), dm_mvet)
dm_mvet[,2] <- c("Active Duty", "Veteran", "Civilian")

dm_soco <- as.matrix(svyby(~ Q11_mr, ~ social, w2_design, svymean, na.rm = TRUE))
dm_soco <- cbind(c(rep(10,2)), dm_soco)
dm_soco[,2] <- c("No", "Yes")

dm_fami <- as.matrix(svyby(~ Q11_mr, ~ family, w2_design, svymean, na.rm = TRUE))
dm_fami <- cbind(c(rep(11,2)), dm_fami)
dm_fami[,2] <- c("No", "Yes")

dm_genr <- as.matrix(svyby(~ Q11_mr, ~ generation, w2_design, svymean, na.rm = TRUE))
dm_genr <- cbind(c(rep(6,5)), dm_genr)
dm_genr[,2] <- c("Silent", "Boomer", "Gen X", "Millennial", "Gen Z")

dm_regn <- as.matrix(svyby(~ Q11_mr, ~ REGION4, w2_design, svymean, na.rm = TRUE))
dm_regn <- cbind(c(rep(5,4)), dm_regn)
dm_regn[,2] <- c("Northeast", "Midwest", "South", "West")

dm_incm <- as.matrix(svyby(~ Q11_mr, ~ income5, w2_design, svymean, na.rm = TRUE))
dm_incm <- cbind(c(rep(2,5)), dm_incm)
dm_incm[,2] <- c("Q5", "Q4", "Q3", "Q2", "Q1")

dm_empl <- as.matrix(svyby(~ Q11_mr, ~ unemployed, w2_design, svymean, na.rm = TRUE))
dm_empl <- cbind(c(rep(3,2)), dm_empl)
dm_empl[,2] <- c("Employed", "Unemployed")

dm_mart <- as.matrix(svyby(~ Q11_mr, ~ married, w2_design, svymean, na.rm = TRUE))
dm_mart <- cbind(c(rep(1,2)), dm_mart)
dm_mart[,2] <- c("Single", "Married")

dm_geog <- as.matrix(svyby(~ Q11_mr, ~ URBAN3, w2_design, svymean, na.rm = TRUE)[2:4,])
dm_geog <- cbind(c(rep(4,3)), dm_geog)
dm_geog[,2] <- c("Urban", "Suburban", "Rural")

dm_educ <- as.matrix(svyby(~ Q11_mr, ~ EDUC5, w2_design, svymean, na.rm = TRUE))
dm_educ <- cbind(c(rep(7,5)), dm_educ)
dm_educ[,2] <- c("Less than HS", "HS", "Some college", "College", "Graduate")

demo_preds1 <- as.data.frame(rbind(dm_mart, dm_incm, dm_empl,
                                   dm_geog, dm_regn, dm_genr,
                                   dm_educ, dm_race, dm_relg,
                                   dm_soco, dm_fami, dm_mvet,
                                   dm_gend, dm_ideo, dm_prty))
demo_preds1$link <- as.numeric(as.character(demo_preds1$Q11_mr)) * 100
demo_preds1$ci_lo <- (demo_preds1$link - 196 * as.numeric(as.character(demo_preds1$se))) 
demo_preds1$ci_hi <- (demo_preds1$link + 196 * as.numeric(as.character(demo_preds1$se)))

fig8_pts <- ggplot(demo_preds1, aes(y = link, x = V1, ymin = ci_lo, ymax = ci_hi, 
                                    label = married)) + 
  geom_linerange(colour = "gray65") +
  geom_point() + 
  geom_text_repel(size = 2.5, direction = "x", vjust = -1.3) +
  scale_x_discrete(breaks = c(1:15),
                   labels = c("Marital Status", "Income", "Employment",
                              "Geography", "Region", "Generation",
                              "Education", "Race", "Religion",
                              "Social Contact", "Family Contact", "Veteran Status",
                              "Gender", "Ideology", "Party ID")) +
  coord_flip() + xlab("") + ylab("Predicted Support for the Military (%)") + 
  ylim(45,85) + theme_bw()
fig8_pts
ggsave("~/Dropbox/Public_Conf_Mil/Chapter 2 Figures/fig2-9_w2.jpeg", plot = fig8_pts, dpi = 700)

#### FIGURE 2.9 - DEMOGRAPHIC PREDICTIONS ####

demo_mod <- svyglm(Q11_mr ~ dem + rep + ideo_m + male + black + hispanic + other_race + income5 +
                     EDUC5 + unemployed + boomer + genx + millen + activeduty + vet +
                     social + family + midwest + south + west + catholic + christian + norelig +
                     city + rural + married,
                   design = w2_design)
summary(demo_mod)

demot <- data.frame(dem = rep(svymean(~ dem, w2_design, na.rm = TRUE), 5), 
                    rep = rep(svymean(~ rep, w2_design, na.rm = TRUE), 5), 
                    ideo_m = rep(svymean(~ ideo_m, w2_design, na.rm = TRUE), 5),
                    male = rep(svymean(~ male, w2_design, na.rm = TRUE), 5), 
                    black = rep(svymean(~ black, w2_design, na.rm = TRUE), 5), 
                    hispanic = rep(svymean(~ hispanic, w2_design, na.rm = TRUE), 5), 
                    other_race = rep(svymean(~ other_race, w2_design, na.rm = TRUE), 5),
                    income5 = c(1:5), 
                    EDUC5 = rep(svymean(~ EDUC5, w2_design, na.rm = TRUE), 5), 
                    unemployed = rep(svymean(~ unemployed, w2_design, na.rm = TRUE), 5),
                    boomer = rep(svymean(~ boomer, w2_design, na.rm = TRUE), 5), 
                    genx = rep(svymean(~ genx, w2_design, na.rm = TRUE), 5), 
                    millen = rep(svymean(~ millen, w2_design, na.rm = TRUE), 5),
                    activeduty = rep(svymean(~ activeduty, w2_design, na.rm = TRUE), 5), 
                    vet = rep(svymean(~ vet, w2_design, na.rm = TRUE), 5), 
                    social = rep(svymean(~ social, w2_design, na.rm = TRUE), 5), 
                    family = rep(svymean(~ family, w2_design, na.rm = TRUE), 5),
                    midwest = rep(svymean(~ midwest, w2_design, na.rm = TRUE), 5), 
                    south = rep(svymean(~ south, w2_design, na.rm = TRUE), 5), 
                    west = rep(svymean(~ west, w2_design, na.rm = TRUE), 5),
                    catholic = rep(svymean(~ catholic, w2_design, na.rm = TRUE), 5), 
                    christian = rep(svymean(~ christian, w2_design, na.rm = TRUE), 5), 
                    norelig = rep(svymean(~ norelig, w2_design, na.rm = TRUE), 5),
                    city = rep(svymean(~ city, w2_design, na.rm = TRUE), 5), 
                    rural = rep(svymean(~ rural, w2_design, na.rm = TRUE), 5), 
                    married = rep(svymean(~ married, w2_design, na.rm = TRUE), 5))
demotpred <- as.data.frame(predict(demo_mod, newdata = demot, se.fit = TRUE))

### Predictions by Demographic Profiles

## Demo 1: Asian female student
demo1 <- data.frame(dem = 1, rep = 0, ideo_m = 2,
                    male = 0, black = 0, hispanic = 0, other_race = 1,
                    income5 = 2, EDUC5 = 3, unemployed = 0,
                    boomer = 0, genx = 0, millen = 1,
                    activeduty = 0, vet = 0, social = 0, family = 0,
                    midwest = 0, south = 0, west = 1,
                    catholic = 0, christian = 0, norelig = 1,
                    city = 0, rural = 0, married = 0)
demo1pred <- as.data.frame(predict(demo_mod, newdata = demo1, se.fit = TRUE))

## Demo 2: Black low-income single female
demo2 <- data.frame(dem = 1, rep = 0, ideo_m = 1,
                    male = 0, black = 1, hispanic = 0, other_race = 0,
                    income5 = 2, EDUC5 = 2, unemployed = 0,
                    boomer = 0, genx = 0, millen = 1,
                    activeduty = 0, vet = 0, social = 0, family = 0,
                    midwest = 0, south = 1, west = 0,
                    catholic = 0, christian = 0, norelig = 1,
                    city = 1, rural = 0, married = 0)
demo2pred <- as.data.frame(predict(demo_mod, newdata = demo2, se.fit = TRUE))

## Demo 3: CA white collar male worker
demo3 <- data.frame(dem = 0, rep = 0, ideo_m = 3,
                    male = 1, black = 0, hispanic = 0, other_race = 1,
                    income5 = 5, EDUC5 = 4, unemployed = 0,
                    boomer = 0, genx = 0, millen = 1,
                    activeduty = 0, vet = 0, social = 0, family = 0,
                    midwest = 0, south = 0, west = 1,
                    catholic = 0, christian = 0, norelig = 1,
                    city = 1, rural = 0, married = 0)
demo3pred <- as.data.frame(predict(demo_mod, newdata = demo3, se.fit = TRUE))

## Demo 4: Hispanic male migrant worker
demo4 <- data.frame(dem = 0, rep = 0, ideo_m = 3,
                    male = 1, black = 0, hispanic = 1, other_race = 0,
                    income5 = 1, EDUC5 = 1, unemployed = 0,
                    boomer = 0, genx = 1, millen = 0,
                    activeduty = 0, vet = 0, social = 0, family = 0,
                    midwest = 0, south = 0, west = 1,
                    catholic = 0, christian = 0, norelig = 0,
                    city = 0, rural = 1, married = 1)
demo4pred <- as.data.frame(predict(demo_mod, newdata = demo4, se.fit = TRUE))

## Demo 5: CA liberal male PhD
demo5 <- data.frame(dem = 1, rep = 0, ideo_m = 1,
                    male = 1, black = 0, hispanic = 0, other_race = 0,
                    income5 = 4, EDUC5 = 4, unemployed = 0,
                    boomer = 1, genx = 0, millen = 0,
                    activeduty = 0, vet = 0, social = 0, family = 0,
                    midwest = 0, south = 0, west = 1,
                    catholic = 0, christian = 0, norelig = 1,
                    city = 1, rural = 0, married = 0)
demo5pred <- as.data.frame(predict(demo_mod, newdata = demo5, se.fit = TRUE))

## Demo 6: NE Boomer woman
demo6 <- data.frame(dem = 1, rep = 0, ideo_m = 1,
                    male = 0, black = 0, hispanic = 0, other_race = 0,
                    income5 = 5, EDUC5 = 3, unemployed = 0,
                    boomer = 1, genx = 0, millen = 0,
                    activeduty = 0, vet = 0, social = 0, family = 0,
                    midwest = 0, south = 0, west = 0,
                    catholic = 0, christian = 0, norelig = 1,
                    city = 0, rural = 0, married = 1)
demo6pred <- as.data.frame(predict(demo_mod, newdata = demo6, se.fit = TRUE))

## Demo 7: Midwestern soccer Mom
demo7 <- data.frame(dem = 0, rep = 0, ideo_m = 4,
                    male = 0, black = 0, hispanic = 0, other_race = 0,
                    income5 = 4, EDUC5 = 3, unemployed = 0,
                    boomer = 0, genx = 1, millen = 0,
                    activeduty = 0, vet = 0, social = 0, family = 1,
                    midwest = 1, south = 0, west = 0,
                    catholic = 1, christian = 0, norelig = 0,
                    city = 0, rural = 0, married = 1)
demo7pred <- as.data.frame(predict(demo_mod, newdata = demo7, se.fit = TRUE))

## Demo 8: Black female NCO
demo8 <- data.frame(dem = 1, rep = 0, ideo_m = 2,
                    male = 0, black = 1, hispanic = 0, other_race = 0,
                    income5 = 3, EDUC5 = 2, unemployed = 0,
                    boomer = 0, genx = 0, millen = 1,
                    activeduty = 1, vet = 1, social = 1, family = 1,
                    midwest = 0, south = 1, west = 0,
                    catholic = 0, christian = 1, norelig = 1,
                    city = 1, rural = 0, married = 0)
demo8pred <- as.data.frame(predict(demo_mod, newdata = demo8, se.fit = TRUE))

## Demo 9: Post-911 white male veteran
demo9 <- data.frame(dem = 0, rep = 1, ideo_m = 5,
                    male = 1, black = 0, hispanic = 0, other_race = 0,
                    income5 = 3, EDUC5 = 3, unemployed = 0,
                    boomer = 0, genx = 0, millen = 1,
                    activeduty = 0, vet = 1, social = 1, family = 0,
                    midwest = 0, south = 0, west = 1,
                    catholic = 1, christian = 0, norelig = 0,
                    city = 0, rural = 1, married = 1)
demo9pred <- as.data.frame(predict(demo_mod, newdata = demo9, se.fit = TRUE))

## Demo 10: Hispanic male marine
demo10 <- data.frame(dem = 0, rep = 0, ideo_m = 4,
                     male = 1, black = 0, hispanic = 1, other_race = 0,
                     income5 = 3, EDUC5 = 3, unemployed = 0,
                     boomer = 0, genx = 0, millen = 1,
                     activeduty = 1, vet = 1, social = 1, family = 1,
                     midwest = 0, south = 0, west = 1,
                     catholic = 1, christian = 0, norelig = 0,
                     city = 1, rural = 0, married = 0)
demo10pred <- as.data.frame(predict(demo_mod, newdata = demo10, se.fit = TRUE))

## Demo 11: NASCAR Dad
demo11 <- data.frame(dem = 0, rep = 1, ideo_m = 6,
                     male = 1, black = 0, hispanic = 0, other_race = 0,
                     income5 = 3, EDUC5 = 3, unemployed = 0,
                     boomer = 0, genx = 1, millen = 0,
                     activeduty = 0, vet = 0, social = 1, family = 1,
                     midwest = 0, south = 1, west = 0,
                     catholic = 0, christian = 0, norelig = 1,
                     city = 0, rural = 1, married = 0)
demo11pred <- as.data.frame(predict(demo_mod, newdata = demo11, se.fit = TRUE))

## Demo 12: NE white collar worker
demo12 <- data.frame(dem = 0, rep = 1, ideo_m = 6,
                     male = 1, black = 0, hispanic = 0, other_race = 0,
                     income5 = 4, EDUC5 = 4, unemployed = 0,
                     boomer = 0, genx = 1, millen = 0,
                     activeduty = 0, vet = 0, social = 1, family = 0,
                     midwest = 0, south = 0, west = 0,
                     catholic = 0, christian = 1, norelig = 0,
                     city = 1, rural = 0, married = 1)
demo12pred <- as.data.frame(predict(demo_mod, newdata = demo12, se.fit = TRUE))

## Demo 13: Male Midwestern farmer
demo13 <- data.frame(dem = 0, rep = 1, ideo_m = 6,
                     male = 1, black = 0, hispanic = 0, other_race = 0,
                     income5 = 2, EDUC5 = 2, unemployed = 0,
                     boomer = 0, genx = 1, millen = 0,
                     activeduty = 0, vet = 1, social = 1, family = 1,
                     midwest = 1, south = 0, west = 0,
                     catholic = 1, christian = 0, norelig = 0,
                     city = 0, rural = 1, married = 1)
demo13pred <- as.data.frame(predict(demo_mod, newdata = demo13, se.fit = TRUE))

## Demo 14: Female Fox News watcher
demo14 <- data.frame(dem = 0, rep = 1, ideo_m = 6,
                     male = 0, black = 0, hispanic = 0, other_race = 0,
                     income5 = 3, EDUC5 = 2, unemployed = 0,
                     boomer = 1, genx = 0, millen = 0,
                     activeduty = 0, vet = 0, social = 0, family = 1,
                     midwest = 0, south = 0, west = 0,
                     catholic = 0, christian = 1, norelig = 0,
                     city = 0, rural = 0, married = 1)
demo14pred <- as.data.frame(predict(demo_mod, newdata = demo14, se.fit = TRUE))

## Demo 15: Male Republican Korean war veteran
demo15 <- data.frame(dem = 0, rep = 1, ideo_m = 7,
                     male = 1, black = 0, hispanic = 0, other_race = 0,
                     income5 = 4, EDUC5 = 3, unemployed = 0,
                     boomer = 1, genx = 0, millen = 0,
                     activeduty = 0, vet = 1, social = 0, family = 0,
                     midwest = 0, south = 0, west = 1,
                     catholic = 0, christian = 1, norelig = 0,
                     city = 0, rural = 0, married = 0)
demo15pred <- as.data.frame(predict(demo_mod, newdata = demo15, se.fit = TRUE))

demo_preds <- rbind(demo2pred, demo1pred, demo3pred, demo5pred, demo4pred,
                    demo6pred, demo7pred, demo8pred, demo11pred, 
                    demo10pred, demo14pred, demo12pred, demo9pred, demo15pred, demo13pred)
demo_preds$ci_lo <- demo_preds$link - 1.96*demo_preds$SE
demo_preds$ci_hi <- demo_preds$link + 1.96*demo_preds$SE
demo_preds$count <- c(1:15)

demo_plot <- ggplot(demo_preds) +
  geom_pointrange(aes(x = count, y = link, ymin = ci_lo, ymax = ci_hi)) +
  coord_flip() + xlab("") + ylab("Predicted Support for the Military") +
  scale_x_continuous(breaks = c(1:15),
                     labels = c("Black low-income single female", "Asian female student", 
                                "CA male white collar worker",
                                "CA liberal male PhD", 
                                "Hispanic male migrant worker", "NE Boomer woman", 
                                "Midwestern soccer mom", 
                                "Black female NCO", 
                                "NASCAR Dad", "Hispanic male marine",
                                "Female Fox News watcher",
                                "NE white collar worker", 
                                "Post-9/11 white male veteran",
                                "Male Republican Korean war veteran",
                                "Male Midwestern farmer")) +
  theme_bw()
demo_plot
ggsave("~/Dropbox/Public_Conf_Mil/Chapter 2 Figures/fig2-9_w2.jpeg", plot = demo_plot, dpi = 1000)

